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Abstract 

We describe a new algorithm for simulating complex Markoff-processes. We have used a 
reaction-cell method in order to simulate arbitrary reactions. It can be used for any kind of 
RDS on arbitrary topologies, including fractal dimensions or configurations not being related 
to any spatial geometry. The events within a single cell are managed by an event handler 
which has been implemented independently of the system studied. The method is exact on 
the Markoff level including the correct treatment of finite numbers of molecules. 

To demonstrate its properties, we apply it on a very simple reaction-diffusion-systems (RDS). 
The chemical equations a+a ^0 and a+b^<6 in 1 to 4 dimensions serve as models for systems 
whose dynamics on an intermediate time scale are governed by fluctuations. We compare 
our results to the analytic approach by the scaling ansatz. The simulations confirm the 
exponents of the a+b system within statistical errors, including the logarithmic corrections in 
the dimension d= 2 . 

The method is capable to simulate the crossover from the reaction to diffusion limited regime, 
which is defined be a crossover time depending on the system size. 


1 Introduction 

Thus the dynamics of chemical reactions in homogeneous systems is well understood, they are 
easy to describe by a van’t Hoff ansatz Q. However, in inhomogeneous systems the transport 
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mechanism has to be taken into account. We want to consider reaction-diffusion-systems 
(RDS) for which the slowest mechanism determines the dynamics of a RDS. RDS are known 
to be capable of a very complex behaviour like formation of spatial and temporal patterns. 
They show this great variety of effects, even if their discrete nature can be neglected and 
they are well described by the mean concentrations n a {r ). Population dynamics, which is 
often described by a formulation equivalent to chemical reactions, deals with much smaller 
concentrations, so the role of fluctuations is much more important than for chemical reactions. 

1.1 The structure of this paper 

At first, to introduce the notation, we briefly review the general method for simulating master- 
equations, with respect to chemical reactions. Then we will present the reaction-cell model 
to describe extensive reaction volumes. After that, we give a detailed description of the 
implementation as a the ” Markoff-automaton”. Our algorithm makes use of non-numerical 
data structures borrowed from computer science, because we need programming methods not 
commonly used in standard numerics. Therefore, we describe our methods in details. They 
may be regarded as ” semi-numerical” algorithms in the sense of ||]. Finally, we will present the 
results we have gained applying our algorithm to the simplest RDS A + A —> 0 and A + B —> 0 
in 1 to 4 dimensions. 

1.2 Reaction-Diffusion—Systems 

We want to denote different types of molecules H a by a Greek index a, while the reactions 
are distinguished by the Latin letter i. The stoichiometric equations describe the chemical 
reaction of the i-th type unambiguously by the initial (or forward) fa and final (or backward) 
b x a coefficients and the related reaction-rate A 1 

S ~ (!) 


The common formulation of the mean-field equations for position-dependent mesoscopic 
concentrations n(x,t) leads to nonlinear coupled partial differential equations 


^-n a = D a A7i a + 

at 


E 


^ - fh) i 


( 2 ) 


In some realisations of certain systems they may also describe the fluctuations @, fjj] or give 
only quantitative changes Ji), nevertheless, it may become necessary to respect the influence 
of fluctuations, especially in biological systems, where the number of individual molecules is 
very small [fj . 

However, in contrast to the simple derivations of the partial differential equations for 
the mean values, the formulation of the equations for the higher moments are much more 
complicated, even for simple systems. 

Note the different scaling of mean concentrations and fluctuations. While concentrations 
are proportional to the number of particles N a in a volume S2, i.e N a = Qn a , fluctuations 
scale with a different power law according to 5N a = \/N a . If we introduce the fluctuations 
of n a as additional quantities 5n a , we have to take into account, that they usually scale 
as Sn a = \/N a . However, the assumption of a fl 1 ' 2 scaling is invalid, if fluctuations may 
increase, as in our case. Other important examples of scaling departing from this law are 
physical systems showing critical behaviour. A rigourous derivation of the fluctuation scaling 
by terms of an 1/fl expansion is given in |y]. 
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Therefore, the analytical treatment turns out to become very difficult, if diffusion and the 
discrete nature of the molecules together must be taken into account. Not only the cost of 
computation grows proportional to the volume, it has to be remarked that it also increases 
exponentially to the statistical moments taken into consideration. 

In the reaction-limited case, diffusion is sufficiently fast to keep the concentrations ho¬ 
mogeneous. Thus the reactions follow a global van’t Hoff dynamics, for which diffusion may 
be neglected. For the diffusion-controlled reactions the time to pass structures like domains 
or empty spaces determines the time scale of the RDS. To introduce spatial resolution we 
subdivide the volume into L d small cells of size ui = Q./L d . 

Therefore, for small numbers of molecules within an uj volume, the discrete nature of the 
molecules must be taken into consideration. The raw estimation for systems of N molecules 
leads to an estimation for the fluctuations of 0{y/N). Thus we need O(N) = 10, 000 molecules 
to guarantee a relative error of percentage order if fluctuations are neglected. 

Former simulations [^| |jj, [uj using cellular automata or lattice gas methods to study the 
dynamics of RDS are restricted to small volumes and small numbers of molecules. Computer 
simulations of complex RDS focus on the mean-field-behaviour taking no notice of the fluctu¬ 
ations, whereas on the other hand stochastical simulations of extensive systems close to equi¬ 
librium are not well suited to dynamics. In contrast, using our algorithm we may handle some 
10' particles and are, furthermore, able to study the crossover from the diffusion-controlled to 
the reaction-controlled limit. Following an analysis of the relevant time scales, the diffusion 
can be simulated on a coarse lattice of cells without affecting the exactness of the results. Our 
goal is to close the gap between solving the mean-field-equations and simulation techniques 
using cellular automata. Both techniques are included, but we want to emphasize, that stan¬ 
dard algorithms may be more efficient to cover these limits. Because our method is related to 
cellular automata, we call it a ” Markoff-automaton” claiming its exactness on the level of the 
Markoff-processes derived. 

The algorithm has been written for systems with a clear separation of the short time-scale 
of reaction events and the long time-scale of the decay of the population. It should be used 
with care, when the scales are mixed, f.e. by local ordering phenomena. 


1.3 The physical systems 

We consider two different RDS consisting of one of the two simplest chemical reactions: 

A + A -> 0, or (3) 

A + B —> 0. (4) 


1.4 The reaction—limited case 


In this case, the diffusion is as fast, that the extinction is controlled by the reaction. Thus 
both systems may be described by a global van’t Hoff ansatz, which leads to the differential 
equations for the concentrations n a 


dt 


UA ( t ) 


d_ 

dt 




The straightforward solution of (H) reads 


-fcni(t), 

— knA(t)n B (t). 


n A (t) 


ua{Q) 
n(0)kt + 1 


(5) 

( 6 ) 

(7) 
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To solve we take advantage of the fact, that the difference A n(t) = riA(t) — ns{t) = An(0) 
is conserved. This way we obtain the mean-held result 


_ nA (0) A n exp(—A nkt) , . 

UA tib (0) — riA (0) exp(—A nkt) ’ 

n B (t) = nA(0) + An. (9) 

For the case of equal initial concentrations, i.e An = 0, equation (^) reduces to (!7;). Note that 
the initial condition ua( 0) = n B (0) is responsible for the algebraic decay. 


1.4.1 The diffusion limited case for A + A —* 0 


In this case the extinction is controlled by the speed of diffusion, which has been analysed in 
detail by the scaling ansatz of Kang and Redner [|j, |hj]. However, at this point, we can only 
give a brief description of their main arguments, because we have outlined their scaling in the 
appendixes. 

While we want to use this simple model to test our algorithm, we also want to examine 
the accuracy of the scaling ansatz, which predicts an algebraic annihilation for the diffusion 
controlled A + A reaction. The time is the time to pass the mean distance between two A 
molecules at t = 0, 


ua(£) oc 


t~ d/2 , d< 2, 
t -1 , d> 2. 


t > t£, 


( 10 ) 


This RDS shows a critical dimension d c = 2, for d> d c the decay is predicted to be determined 
by mean held behaviour. 


1.4.2 The diffusion limited case for A + B —> 0 


In the A + B reaction we obtain the time scale determined by the time to pass a domain as 


1 

k= D 


\/ua( 0) - \Jn B (0) 


-2/d 


( 11 ) 


An algebraic annihilation uab oc t * with d c = 4 appears on a time scale t < t%. Thus, if 
iia( 0) < n B (0) we obtain in the limit An = nA(0) — n B (0) —> 0, 


n,A, B (t ) oc 


t~ d/i , d < 4, 
£ -1 , d > 4. 


t < tt, 


( 12 ) 


A more detailed explanation of the scaling arguments is left to the appendixes. Both of these 
simple reactions are the most important examples for the influence of fluctuations in RDS. 
Therefore, to improve our understanding of RDS beyond mean-held-behaviour, we present 
our simulation method as an approach correctly treating dynamics and statistics in a natural 
way. 


2 Simulating master-equations 

We want to treat a chemical reaction as a Markoff-process, which is described by the time 
evolution of the probability Pjt to be in state X at time t following a stationary master- 
equation 


d Px 

dt 

= £ 

Px- 

W x^X' Px} 


X' 

inhow into X 

outflow from X 
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in our case the discrete vector X denoting the number of molecules |I^ ]. 

The minimal-process-method simulates the master-equation as a random-walk in the 
space of all possible configurations. Starting the random-walk in state X, we only need to 
determine the time leaving X, and the successor state X ', which is equivalent to draw the 
lifetime r g and the transition X —> X' as random numbers. 

For both steps the knowledge of the flux out of A' is completely sufficient. The lifetime 
rg is exponentially distributed, while the selection of the transition requires the drawing of 
a random-number X' proportional to the transition rate Wx^x' ■ This property is due to 
discrete Markoff-processes, and we want to emphasize that it is exact and not affected by any 
further assumptions. 

To explain the name ’’minimal-process-method” it has to be remarked, that the master- 
equation (^d|) may be interpreted in different ways. We may look upon the outflow term as the 
death-process and upon the flux into a state as the birth-process of the state X □•0- The 
minimal-process being the process with the least number of events is the random-walk -process 
in the discrete configuration space of all possible X. The minimal process distinguishes on the 
rule, that a death process in X always coincides to a simultaneous birth process in X'. With 
the sum 

(14) 

x' 

the way to generate this random-walk and the sequence of states X ( t ) is given by the following 
simple algorithm: 

Minimal—process—method 

1. Generate an exponentially distributed random number r with 
Pr(r) = Wg ■ exp(— Wg ■ r), let t <— t + t, 

2. select X' according to the probability 
Pr(A^A ') = W^,/W^ 

3. go back to step 1. 

The characteristics of discrete processes were already known to Markoff, the minimal-process- 
method at least since the early sixties |g|, . However, no one has made use of it for computer 

simulations until the mid seventies when Gillespie applied it to chemical reaction-systems 

6S- 

2.1 Chemical reactions in a homogeneous volume 

For the present, we want to discuss arbitrary chemical reactions in a small volume u>, for which 
the diffusion is so fast that any spatial inhomogeneity may be neglected. After that, we will 
introduce diffusion as a reaction like process of molecules leaving the volume u>. In this way, 
we want to describe the dynamics by which diffusion relates distances of an extensive length 
scale to a time scale. 

According to (jlj) the number X a of molecules E a changes from X a to X a +b l a — each time 
the reaction of the i-th type occurs. I.e. /), molecules E a are consumed and b 2 a molecules E a 
are produced. Forward and backward reactions shall be distinguished by their enumerations. 
If the reactions j and i are reciprocals of each other, the role of the coefficients is exchanged, 
i. e. f a = V a , and f 3 a = b l a . 

Since the number of molecules is an extensive quantity, the number of chemical reactions 
per unit time must be extensive, too. The rate of events, the reactivity A* of a chemical 
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reaction is defined as 


j ("the reaction i occurs duringi 

^ r (the time interval [f, t + dt\ J ' 

According to van’t Hoff [Q], it is given by 


XJ 


{X a -/*,)! w/i 


(15) 


(16) 


where A* is an intensive constant which does not depend on X a . Note that the reactivities A 
depend on the state, i.e. A 1 = A*j, A = A-^, the index X is usually suppressed. The fraction 

XJ = X a (X a — !)•■• (X„ -f*+l) 

{X a -fi)lufi u& 

takes into consideration that each individual molecule can only be consumed once. This is 
important, if the numbers X a are small. For large X a , this fraction turns into the more 

familiar form ignoring any power up to *. 

The van’t Hoff approach suggests chemical reactions being Markoff-processes [Tl[ 
which may be justihed due to a stosszahlansatz for homogeneous concentrations. The config¬ 
uration is the vector consisting of the number of molecules 

X = (AT, A' 2 ,..., X a ,...). (18) 


For abbreviation we group the initial and final coefficients for the reaction of the i-th type as 
vectors /* = (/(, f %,..., /J...), b l = (b\, b l 2 ,..., b l a ,...). Each type of reaction is realized 
by an integer vector operation, the vector 8 l = b 1 — /* denoting the changes caused by the 
i-th reaction: 


X <- x + s\ 

(19) 

In our case the vector X is simply a one- or two-dimensional vector 


- f (Aa) reaction A + A — > 0, 

( Xa,Xb ) reaction A + B — > 0. 

( 20 ) 

The reaction-rates are computed according (|I(j) as 


. , Aa(Aa - 1) 

AAA(aJ » 5 

UJ Z 

( 21 ) 

. , X A X B 

AaB = a ABU 

UJ Z 

( 22 ) 

The homogeneity of the concentrations within u is necessary for the van 
will have to subdivide the volume into sufficiently small cells to justify a 
introducing diffusion as hopping from one cell to an adjacent one. 

’t Hoff ansatz. We 
local form of (|l 6 |) , 


2.2 Gillespie’s algorithm 

With these assumptions Gillespie describes the sequence of states of a chemical system by the 
random-walk of the vector X = X(t). His algorithm is a version of the minimal-process- 
method adapted to the dynamics of chemical reactions, which is based on the well known 
properties of discrete stochastic processes|hi|, [IT[ [Tf]. |7^|. We summarise his main results 
according to our requirements. The total reactivity is the rate 

A = 51 A "’ ( 23 ) 

all reactions i 
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which determines the probability of a chemical reaction during the infinitesimal time interval 
[t, t + dt] according to 


A dt = Pr 


{ some chemical reaction changes the 
state X during [t,t + dt] 


(24) 


A is the inverse of the average lifetime of X, therefore the random time ta until the next event 
occuring has the probability density 


Pr{rA,time to next reaction } = Aexp(— Ata). 


(25) 


Because we have assumed that chemical reactions are Markoff-processes, this equation yields 
independently of the time when the last change occurred, which is a general statement for any 
discrete Markoff-process. According to Gillespie, the probability of the reaction of kind i is 
proportional to its contribution to A, thus 

A* 

Pr{Reaction i} = — (26) 

is the probability, that the next change of X in the chemical system is determined by a reaction 
of type i defined by ©■ The transition rate is given by the probability of all chemical 

reactions leading from X to X' 

= E Ai - ( 27 ) 

all reactions i with 
X' = X — f + hi 

Note that different chemical reactions may change X by the same <5 l = /* — b l because only 
the backward and forward stoichiometric coefficients f l and b 1 together determine the reaction 
unambiguously. 

Gillespie’s algorithm implicitly uses equation ([ 27 ]) to decompose the relation of reaction and 
transition probabilities. For arbitrary chemical reactions it computes the random time-step ta 
until the next reaction by drawing an uniformly distributed random number rnd € [0,1). The 
reaction type i usually is selected by a simple loop, known as linear selection algorithm using 
the sum s for the integration. The variable t denotes the time, proceeding in exponential time 
steps. The order of reactions denoted by i is arbitrary. After the reaction has been carried 
out, the reaction reactivities A 1 and the total rate A = W^ are computed again. 


Gillespie’s algorithm 

1 . ta <- log(l — rnd), t <— t + ta 

2. r <— rnd, s <— 0, i' <—first reaction, 


3. while s < rA, 

(a) 

* <— i , i' <—next reaction 


(b) 

add up s <— s T A* 


4. do reaction type i by X <— X + <P, compute A 1 , and A <— A 1 , 

5. go back to step 1. 


Note the adaptation to the time scale A -1 , which is responsible for the high flexibility of 
Gillespie’s algorithm. Thus it is a fast and easy to use Monte-Carlo method for the dynamics 
of arbitrary reaction-systems and has the main advantage over the mean-field description, 
that it takes into account internal fluctuations caused by the finite numbers of molecules. 
Because the random-walk imitates the sequence of states of the simulated Markoff-process, 
any correlations and higher moments may be measured as in a real experiment. 
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3 The reaction—cell method 


The molecules are treated as point-like particles with a finite interaction probability, impli- 
cicelty assuming, that their Brownian motion is independent. Our simulation introduces a 
reaction-cell method. In general, for small diffusion-constants, the homogeneity condition of 
the van’t Hoff ansatz cannot be fulfilled for the total d-dimensional volume 57. Thus 52 is sub¬ 
divided into small cells (cubes or squares etc.) of volume u> — ■£[ = h d . We have to impose the 
condition, that the size of uj may be justified by a local van’t Hoff ansatz. This prerequisite 
depends on the time-scale tr between two subsequent chemical reactions compared to the 
time scale tb of a molecule leaving a > by diffusion. The diffusion-time-scale has to be much 
faster than the reaction time-scale, i. e. tb -C tr with td = h 2 / D and tr = 1/A. 

Our algorithm covers the held between the deterministic description of RDS, which does 
not describe the fluctuations correctly, and the cellular automaton approach, which is not well 
conditioned to describe large numbers of particles. The assumptions have to be justified in any 
case, and have to be modified carefully, f.e. if the spatial extension of the moving objects has 
to be taken into account. This is important for the A + A —> 0 reaction in lower dimensions 
00- where the seed reaction constant may become a irrelevant quantity. In this case first 
principle Monte-Carlo simulations have to be performed and our algorithm has to be checked. 

3.1 Reaction-Diffusion—Systems 

We distinguish different cells by an integer vector index n. Therefore, the possible reactions 
now depends on the local situation which is described by the number X af t of molecules H a in 
cell n. The total volume 52 is represented by the direct sum 

^ = ©^ s = ®®-^“ s ’ ( 28 ) 
n fi ex. 


consisting of sub-vectors Xft = {. .., X af i, ...} of the numbers of molecules in each cell. 

For the A + A —> 0 system X has the form 

X = { (JCa (1 . 1} ) , (Xa ( i,..., 2) ) , • ■ ■, (AA (i ,...,i)) } , (29) 

respectively for the A + B —> 0 reaction 

X = {(AOua,... u), -Xb(i.,., A b(1i .,., 2) ) , (30) 

■■■, (Xa(l,...,l'),Xb(l,...,l))} 

In this section we want to introduce diffusion as a random-walk on the lattice of reaction- 
cells. The diffusion in an arbitrary, not necessarily euclidian, topology may be regarded as a 
reaction-like step of a molecule H Q from n to one of its next neighbours m, 


For a diffusion-step of a molecule a the rate 


= D a 


1 

! 7T 


with .Dafjm denoting the local diffusion-constant, describes the probability 

f a single molecule E a jumps from n to) 


A 0 


\dt = Pr • 


m during the time intervals [t,t + dt] 




(31) 


(32) 


(33) 



On a lattice built of volumes uj = h d , we obtain in the limit t —> oo the correct Brownian 
motion behaviour by the central-limit-theorem. The symbol (n) denotes the set of all 2 d next 
neighbour cells of n. Due to the symmetry of the rate to do a single step into any direction 
i.e. the rate of leaving n by diffusion, is D af t/h 2 , with 


A aft — 




Acmm 2dA 0 


(34) 


In the isotropic and homogeneous case we omit the spatial indices Ac = X a n and the rate for 
stepping in a certain dimension gets A a /2 d, while the total rate for leaving a cell is \ a . 

If we have X ai n molecules H Q in cell n, the rate for the next diffusion step of any E a to rh 
is 

Aftrn — ^nrnXfxn. (35) 

This rate is the equivalent formulation to ©> treating a diffusion-step like a first order 
chemical reaction. If a diffusion-step occurs, another cell rh gets involved, because one E Q 
moves from n to rh. Thus a diffusion-step is carried out by 


1 . 


2. X, 


X a rh T 1. 


(36) 


Notice that in our case A a Hfh depends only on the contents of n, from which a molecule 
hops into a neighbour rh £ (n), although the assigned diffusion-step changes the contents of 
the other cell rh, too. Thus a diffusion-step changes the state of two cells, however, its rate 
depends only on the contents of the first cell n triggering the event. 

For the chemical reaction of type i we simply have to change the contents of h by Xn <— 
Xft + 8 l . Therefore, the probability Qn of an event triggered by cell h is the sum of reaction- 
rates Ax and diffusion-rates A a nrh 

Qn=J2A i n + Y^ A “™ ( 37 ) 

i ol rh£ (fi) 

Because all reaction- and diffusion-steps are assigned to the cell n, we may identify an arbitrary 
step entirely by an integer i and the cell n. We may ignore the additional index rh for the 
selection mechanism. A reaction may be defined by an addition of its associated difference 
<5k, i.e. X <— X + 5^, involving one or two cells. In section 6.1.1 below we suggest obvious 
generalisations of rates depending on the contents of the neighbourhood of h. This way it is 
easy to introduce dynamics depending on gradients, like the movement of kinks and steps on 
surfaces. 


3.2 Times scales of reaction— and diffusion—processes 

If we want to fulfill the assumptions of a local van’t Hoff ansatz we have to satisfy the condition, 
that diffusion is much faster than any chemical reaction. More precisely spoken, this means, 
that the local reactions A)j and diffusion- A Q ^ rates define different time scales, which have 
to be separated. Therefore we choose the size of our cells sufficiently small, for the condition 

An < A aft rn, for any i, a (38) 

to be satisfied. This can always be achieved, because reaction probabilities are extensive 
quantities Ab oc uj = h d , while diffusion-rates increase with h ~ 2 . Thus A a n.rh oc u}~ 2 ^ d and 
Aanm oc uJ 1 ~ 2 ^ d , the lattice constant h always can be reduced to make the ratio 

-r^_ oc u} 2/d = h 2 (39) 

cxfirh 

arbitrarily small. For any practical purpose it is not possible to give an a priori estimation, so 
we recommend to choose uj as large as possible rejecting the simulation runs if 10Ak > AL-. 
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4 Selecting a single cell by the method of logarithmic 
classes 


Although we have argued that there is no principal difference between the Gillespie algorithm in 
a homogeneous volume and the method on a lattice of reaction cells, a significant complication 
arises as the number of possible changes in each Markoff-step is an extensive quantity. If we 
make a straightforward generalisation of Gillespie’s algorithm using a linear selection strategy 
we run into the problem of selecting a single event from an extensive quantity of transitions. 
The total reactivity 

q=E Qb = EE q « (40) 

ft ft i 

again determines the mean time step. According to condition ( ps| ) we choose the sequence of 
transitions considering at first the more probable diffusion-steps, afterwards the possibility of 
chemical reactions. 


Gillespie’s algorithm for reaction—cells with linear selection 

1 ■ tq < -ilog(l-rnd), t^t + r Q , 

2 . r <— rnd, s <— 0, 

3. while s < rQ, for all cells n , 

(a) next n, 

(b) while s < rQ, 


i. next i, 

ii. sum s - 


s+Ql. 


4. do transition X t— X' + 8~, 

5. for all cells n and reactions i involved, 
recompute all Q' n => Qft => Q, 


6 . go back to step 1. 

This method suffers from the severe disadvantage that the sum S computed at step (|3) 
adds up an extensive number of reaction probabilities. Drawing a random cell due to its 
contribution Qft to Q, the loop consumes computing time proportional to the number of cells. 
This procedure is completely unsuitable for a single step of a computer algorithm, see figure 


If we do not have any additional information, we have on the average to walk the half 
length of the main loop. Even for a small lattice of 10 4 = 100 x 100 cells this means an 
increase for the computing time compared to the algorithm without spatial resolution by an 
average factor of 5, 000. It seems to be impossible to circumvent this problem by rearranging 
the cells or some kind of pre-sorting, because the selected reaction changes the cell reactivity 
every time. A random selection according A. J. Walker’s Q algorithm would be efficient only 
if the reaction-rates did not depend on time. The extensive quantity of cells to be considered 
is a principle obstacle. 

These reflections show that linear selection is totally inefficient. Thus we have tried another 
algorithm based on the von Neumann rejection, which selects at first the reaction-cell and in 
a second step the reaction or diffusion in this cell. 
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Figure 1: Selection algorithms: (a) Linear s., the loop integrates Qn up to rQ , thus the cost is 
extensive, (b) Von Neumann rejection without improvement, the cost is the ratio of the area of 
the rectangle marked off by the upper limit Q to the area below Q A single peak may severely 
affect its efficiency, (c) s. by the method of logarithmic classes, the upper and lower limits for each 
class are denoted by dotted lines. The acceptance a > 0.5 is the ratio of the reordered Qn to its 
upper estimation 2^ 001 (Id(Qjj)) + 1, we find the cost being approximately constant. 
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4.1 Von Neumann rejection 

The von Neumann rejection requires an upper bound for the probability to select a cell. 
Therefore we have to make the assumption, that 


Qft < Q for each cell n. 


(41) 


We replace step (|s|) by 


§ 

von Neumann rejection 


(a) select a cell n with uniform probability, 

(b) choose a uniformly distributed random number r <— rnd, 

(c) if Qn > rQ then select a reaction according to X af t 

else reject this choice and go back to step (a). 





On the first glance this algorithm seems to solve all problems. The method does not depend 
on the number of cells, thus the problem of the dependency on the extensive number of cells 
does not arise. However, a more precise look reveals that the efficiency of the von Neumann 
rejection strongly depends on the homogeneity of the RDS. This is demonstrated in figure 
□(b). The acceptance ratio a is the quotient of the area below Qa and the rectangle delimited 


by Q 


a = 


XX Qft 

QZn 1 ’ 


(42) 


and 1/a is proportional to the number of runs through the loop needed to select a cell. 
Therefore the cost C v nr of the algorithm is proportional to the inverse of the acceptance 
ratio, 


I QV-1 

WNE (X — — -=- 

II XX *3™ 


(43) 


For very inhomogeneous systems, i.e. for systems dominated by a single peak, this method 
may slow down by an arbitrary factor, which has been studied for the simulation of a biological 
system [[|. This disaster can be avoided, if we split up the cells into classes according to their 
dual order of magnitude. The acceptance ratio is improved to 75%. 


4.2 Method of logarithmic classes 

If we want to handle all orders of magnitude of the cell reactivities Qn we have to implement 
a logarithmic classification scheme. We define the logarithmic class L z as the set of all cells 
n with a reactivity Qft in the same order of dual magnitude. The symbol ld(a;) = log 2 (:r) 
denotes the logarithm with base 2, and the function floor(a;) denotes the largest integer which 
is not greater than x 

L z = {n\ floor(ld(Qjs)) = z} ■ (44) 

The reactivity 4 of a certain class is given by the sum of the reactivities of its elements 

4 = Q*- ( 45 ) 

n£L z 

There is no principal restriction of the range of 2 

2 = floor(ld(Qn)) 6 {-oo,..., -2, -1, 0,1, 2, 3,...}. (46) 


12 









The class L_oo contains all cells without any possibility of a reaction- or diffusion-step, above 
all the empty cells. Because L-oo cannot be selected, i.e £-00 = 0, there is no need to 
represent it in computer memory. For any practical application, 2 is restricted to a finite 
range 2 m i n < 2 < 2 max . The intersection of two different classes is empty L z fl L z i = 0, 
z ^ z' ', thus the total reactivity is expressed by the sum 


00 2 max 

Q = ^Qn = Yh iz= Z *'■ 

n z=-oo z=2 min 


(47) 


The inequality 

2 z <Q h < 2 2+1 , He L z (48) 

guarantees, that the reaction-rates are relatively homogeneous within a certain class. The 
rates of two members of any class do not differ by more than a factor of 2 


max Qfi 

nEL z 

min Qft> 

n'eL z 


< 2 . 


(49) 


This is an ideal starting point for the von Neumann rejection. We therefore decided to imple¬ 
ment an event handler, which is able to select a reaction-cell according to its reaction-rate by 
a von Neumann rejection step within the class. The number of classes is small, i.e. of 0(10), 
thus the choice of a class is based on a linear selection algorithm. 


4.2.1 Selecting a reaction 


The problem has been split into three qualitatively different steps: 

1. select a class L z with probability £ Z /Q, 

2. select a cell n £ L z with probability Qn/£ z , 

3. select a reaction-step i within cell n with probability Q^/Qn- 

We stress again, that the index i represents both reaction and diffusion transitions. The last 
two steps are selected according their conditional probabilities 


Pr{select class 2 } 
Pr{select cell n | class 2 has been selected} 
Pr{select reaction i | cell n has been selected} 
therefore the probability of selecting a reaction in a cell 
Pr{select reaction i in cell n} = 


4 

Q’ 

_ Qfi 

TI' 

Qk 

Qn 


On 

Q 


(50) 

(51) 

(52) 


(53) 


has been maintained correctly. 

The algorithm choosing a class is a linear selection. However, it is not necessary to order the 
sequence of classes ( Zi ,..., Zj) in a naive way with z z < Zi+\ or 2 * > 2 j +1 for all i. Moreover, it 
has a favourable effect, if the classes are sorted with respect to their reactivity, i.e £ Zi > 4 i+1 - 
To speed up the selection, the classes of the highest probabilities to be selected are successively 
moved to the head of the loop by a bubble-sort, see step ( pc[ ) below. This is the most efficient 
algorithm to select a reaction in a class according to its rate £ Zi , because the classes with the 
smallest probabilities to be drawn least are moved to the tail. 
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1. Selection of a logarithmic class L z 

1. random number r <— rnd, sum s <— 0, index i = 0 

2. while s < rQ 

(a) increase i <— i + 1, 

(b) add up s <— s + l Zi , 

(c) if i > 1 and £ Zi > £z i _ 1 then 

exchange the order of the classes Zi and z z -i. 

3. return z = Zi, i.e draw L Zi . 

After z has been selected, it is easy to propose a cell n by drawing a uniform random 
number u £ {0 ,... ,v z — 1}. We have implemented each class L z as an array F z [0,..., v z — 1] 
of v z =| L z | elements, each describing the state of a single cell. The following subroutine does 
the von Neumann rejection of a cell in the previous selected class z. 

2. Selecting a cell n in class L z by von Neumann rejection 

1. propose u <— floor(i/ 2 rnd), n F z [u], 

2. draw a uniform random number r *— rnd, 

3. if Qn > r2 2+1 then return n, 

else go back to step 1. 

Because of the inequality (^) we know lower and upper limits of the reactivity, 2 Z < Qft < 2 Z+1 . 
Therefore, it is guaranteed, that the probability for a rejection is less than 0.5. If we assume 
that the rates in each class are distributed uniformly, we get an acceptance ratio a = 0.75. 
In the cell n the reaction is chosen by a linear selection method, because the number of 
possibilities usually is small 0(1) ... 0(10). 

3. Selecting a transition i within cell n, 

1. random number r <— rnd, sum s <— 0, 

2. for all diffusion-transitions j: 

(a) sum s <— s + Ql, 

(b) if s > rQn then return step j, 

else next j, 

3. for all chemical reactions i, 

(a) sum s *— s T Qh, 

(b) if s > rQn then return step i, 

else next i 

The changes in all cells and classes involved have to be registered by bookkeeping steps which 
require the computation of the reaction-rates Qh, the local reactivity of a cell Qft, its class 
reactivity £g and the global reactivity Q. For a diffusion-step concerning a further cell m the 
bookkeeping has to be done for these cells, too. 
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4. Bookkeeping for all cells n, m involved 

For all cells n, as far as their rates are concerned do 

1. compute the local rates Qn, 

2. compute the new cell reactivity Q'n and its new class 
z’ <- floor(ld(Q' i )), 

3. if the class has not changed, i.e. z = z', 

then update i z <— l z — Qn + Q’n-> 

else delete cell n in class z and insert it into z ', 

including an update i z <— £ z — Qn and l z > <— l z , + Q'n- 

4. update the total reactivity Q <— Q — Qn + Qn and the cell reactivity Qn <— Q'n> 


At this position the algorithm has been presented except for the data structure. Inserting and 
deleting of cells is somewhat delicate and will be the subject of the implementation section. 

4.2.2 Implementation of the event handler 

Our approach makes use of a lot of sophisticated, non-standard and non-numeric algorithmic 
structures. Because we do not see any way to obtain these results by usual programming 
methods, we want to introduce our data structures in details. The procedures dealing with 
these structures have been presented in the previous section. 

Our program has been developed in C to achieve the highest portability and speed. We 
have written the code in an object oriented style nevertheless using only the ANSI C 2.0 
standard. We implemented the following duties 

1. an initialisation procedure, which serves as a constructor, to define the numbers of events 
and classes to be managed, 

2. a procedure, which can be seen as a destructor, to clean up the data structures, 

3. a random generator drawing a cell according its reactivity due to its contribution follow¬ 
ing the algorithm of the section above, 

4. a subroutine to insert an event, 

5. another subroutine for the deletion of an event. 

4.2.3 Basic data structures 

The most flexible data structures are pointers, thus an event is described by a pointer to the 
cell being represented. Because we do not want to make any assumption about the cell itself, 
we have implemented an event as pointer to the empty type, i.e. the type event is equivalent 
to the type (void *). This choice has the advantage that it may easily be converted to a 
more problem-adapted structure. The event selection has been kept completely separated to 
encapsulate the problem-independent structure of the event handler from the problem-related 
structure of the cells and the topology. 

The cell-structure contains all informations about the contents of the cell, in our case the 
integer numbers of A or A and B molecules. Unfortunately we found that the number of events 
contained in each class changes rapidly in unpredictable way. It is not possible to specify an 
a priori strategy to determine the size of a certain class in advance. We have therefore been 
forced to write an reorganisation procedure balancing the memory by dividing it into parts 
proportional to the size of each class. 

The link event —> cell is used for the selection mechanism, i.e for events caused by a 
reaction of n. This may be a chemical reaction within n or a diffusion-step originating in n. 


15 






Figure 2: The cell-event data structure 


The link from cell —> event is necessary to respond to changes in cell ft triggered by other 
cells to concerning ft, i.e. some molecule diffusing into ft. 

Since the event-structure may be moved by memory reorganisations, the cell —> event 
link has to be updated every time an event is moved. This is achieved by addressing via the 
link toward the other direction event. cell. e <— event. In general, the structure moving has 
to inform its related partner about its new address. 

The cells and the events are stored in two separate arrays. The cells are arranged in a one 
dimensional array of size L d using helical boundary conditions. The topological neighbourhood 
of a single cell belonging to a certain index array n= ^!s=o L s ns with 0 < n < L d — 1 can be 
addressed by 

n ± L s modulo L d , 0 < 8 < d — 1. (54) 

This way the topology has been implemented independent of the integer dimension d. For 
chemical reactions in fractal dimensions a more sophisticated structure is required. In this 
case the neighbourhood of a cell has a much more complicated topology. 

The memory overhead caused by the double link between event and cell data structures is 
small and becomes less important with an increasing number of reacting species j7j. 

4.2.4 Data structure of the logarithmic classes 

The array of events has a substructure appropriate to the contents of the logarithmic classes. 
For each class a descriptor contains all related information. 

class{ 

pointer to next class L Zi+1 ; 
integer z z , 
derived information 
double float l z -, 
integer v z -, 

F 2 [] as a pointer to e ; 

} 

Within each class the events are stored in a one dimensional array, which is to contain no empty 
places. The classes are realized as structures containing several elements, as shown above. The 


sequence of classes as a linked list 

the power of 2 related to the class 

for speedup only 

class reactivity 

number of events in class 

the first element of the class 
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sequence for the linear selection is realized by a linked list, so there is a pointer to the next 
class. The integer z contains the power of two managed by the structure. For speedup, some 
related frequently used information as 2 Z ,2 Z+1 etc. also may be stored in this structure. The 
number v z describes the variable size of the array of events F z [], which is implemented as a 
pointer to its first element. At this point, we use the ability of the programming language C 
to access entries by indirect or relative addressing. The complete data structure including the 
links to the cells is shown by figure |^. 

If an event is inserted, it is simply stored at the end of the (sub) array F z [], after that 
v z <— v z + 1 is increased. If an event i is deleted, gaps in the array of events F z [] must be 
avoided to keep the von Neumann rejection efficient. Therefore the last event of the sub-array 
has to be moved in the gap, F z [i] <— F z [u z ], before decreasing v z <— v z — 1. Because there is 
a cell structure pointing to this event, its cell has to be informed that its related event has 
moved. 

Since there is no a-priori strategy to estimate the strong fluctuations in the size of the 
classes, we run into the problems of collisions in memory. If the array of events F 2 [] of a 
single class L z is going to overwrite the array of events F z +i[] of its successor, it is necessary 
to reorganise the division of the memory occupied by the array of events. This is done by 
crunching and re-expanding all classes proportional to their actual memory consumption. 


4.3 Implementation of the reactions studied 


According to our approach to separate the event handling from the chemistry within a cell, 
most of the work is done once the event ha ndler h as been implemented. The cells are realized as 
structured variables according to section (4.2.4). Thus adapting the algorithm to any specific 
RDS does not require more than the implementation for drawing a reaction within a cell to 
select the specific reaction. In our case we have to choose 


1. a volume Q, 

2. the number of cells L d defining a discretisation uj = Q/L d and h = a j 1 ^, 

3. a hopping rate \ a ftrh = D a /h 2 , 

4. and a local reaction-constant A*, = uik/ui 2 = k/uo. 

Another procedure has to be provided to do the related changes and to compute the reaction- 
rates. According to (jbij) and (|35|) the reaction- and diffusion-rates of cell n are given by 


An 

= XkX A n {XaH 1) 

or 

(55) 

A fi 

A/c X A nX Bfi 

and 

(56) 

A a ft 

XoL'fiXot.rirh 

with a = A, B. 

(57) 


5 Simulations in several dimensions 

The aim of our simulations is to test the accuracy of the scaling ansatz. We want to examine 
both limits of the reaction-diffusion-system. On the one hand, there is the reaction-controlled 
limit, where the dynamics is determined by the reaction-constant D/D 2 ' d k , expecting 
mean-field behaviour n a (t) oc t -1 , on the other hand we examine diffusion-controlled limit 
D/D 2 ' d -C k, where the long time annihilation is delayed by the necessity of transport. 

Defining length and time scales we have chosen S2 = 1 and D = Da = Db = 1. This way 
the reaction constant k is fixed, too. Since k now has the meaning of a rate in the case of one 
molecule per unit volume, D and k get inverse times [ D\ = [fc] = [1/t], 
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classes events cells 

Figure 3: Data structure of classes, events and cells. At first a class is selected, then an event 
related to a cell, and finally a reaction in a cell. The class structure points to the first of its events. 
Among the event sub-arrays of different classes empty places are left. 
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Figure 4: Formation of domains in the diffusion-controlled reaction A + B —» 0 for a volume of 
100 x 100 cells. The symbol ” I” denotes cells containing one or more A molecules while the symbol 
denotes B. 
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Figure 5: Simulation of A+B —> 0 for diffusion-controlled reactions in the dimensions d = 1,2, 3,4, 
extension L d = 3 ■ 10 4 , 400 2 , 50 3 , 25 4 , An = 0 


5.1 The system A + B —> 0 

Figure (Q) shows a snapshot of the formation of domains, figure (Jo]) shows the time dependent 
density n a ( t ) for the diffusion-controlled A + B —> 0 reaction. The logarithmic plot compares 
the results of our simulations to straight lines t~ d A predicted by theory. The size of the 
system obviously is sufficiently large, thus one run in each dimension may be considered as 
self-averaging. The figures (pj) and (|7j) show the exponents of the algebraic annihilation for 
both diffusion- ( D <C k) and reaction- ( D k) controlled limits. In our simulations these 
inequalities are realized by a ratio of at least 9 orders of magnitude for k D = 1. To be 
sure to be sufficiently close to the limit k/D —> oo the rate k has been increased until the 
exponents measured by linear regression have not changed any more. 

Obviously, the scaling theory for the A + B —> 0 reaction is within the statistical errors of 
a few percent. 

5.2 The system A + A —■» 0 

For the A + A —> 0 reaction we observe logarithmic behaviour in its critical dimension, as 
predicted by |l^|. This is not a multiparticle effect and results from two dimensional Smolu- 
chowsky theory for only two particles The has divided by log(t) to reproduce the correct 

t~ l scaling. The results have been obtained by k = 10 s = 10 9 ... 10 12 , thus we are convinced, 
that our simulations give an excellent approximation to the diffusion-controlled limit. The 
simulations have been done with some 10 7 molecules. For every case studied, we have per¬ 
formed a finite size controi by doubling the size L until we have not find any influence of the 
finite volume. Furthermore, we could neglect corrections, even if in d = 4 where we have been 
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A + B —* 0, exponents 7 
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Figure 6: Simulation of A+B —> 0 for diffusion-controlled reactions in the dimensions d = 1, 2, 3,4, 
numerical values in table (Q) 
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limited to a hypercube of linear extension L = 25, what may be a consequence of d c = 4. 


6 Conclusion 

Thus within statistical errors of 1. .. 2% and an uncertainty due to the method of linear 
regression of the same order of magnitude, the tables (Jl[) and (^|) show, that scaling theory 
describes the A + A —> 0 and A + B —> 0 correctly, including the important cases A+A —> 0 and 
the logarithmic corrections in d = 2, where we found an annihilation according to riA(t) oc 

Thus, we have shown that for the simplest reaction-diffusion-systems, in which the influ¬ 
ence of the finite number of molecules in a single cell must be taken into account, the algorithm 
is able to cover the whole field from diffusion- to reaction-controlled systems. A detailed study 
of the crossover of both regimes which determines the crossover time depending on the size of 
the system and the diffusion constants will follow. For the simplest reactions, the theoretical 
results have been reproduced nearly exactly. 

Our approach follows a strategy opposite to the common cellular automaton philosophy. 
Instead of simulating a lot of simple diffusion-steps represented by integer operations, we try 
to find a length scale below which diffusion may be neglected. On this scale one single hopping 
process may replace a large number of infinitesimal automaton lattice gas steps. 

The price to pay is the implementation of a complex event handler to maintain the exact 
Markoff-probabilities. From our experience we are pleased that the total overhead produced by 
the event handler is less than a factor of 3-4 per Markoff-event compared to reaction-systems 
without diffusion. The requirement of diffusion being much faster than chemical reactions 
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Figure 7: Exponents of the algebraic decay of A + A —► 0 for different dimensions and several 
reaction-probabilities k = 10 s , numerical values in table (j|) 
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d 

1 

2 

3 

4 


0.250a 

0.534a 

0.8264 

0.9786 

Is 

0.256a 

0.512a 

0.803s 

0.9997 




0.799s 





0.7826 



0.249i2 

0.532n 

0.768a 

1.013 9 

7oo 

0.243i2 

0.540n 

0.759a 



0.246i2 





Table 1: Numerical results of A + B —> 0, see figure (j|). The values are presented as j s , where 
the subscript s denotes the order of magnitude of k = 10 s with 700 representing the best approach 
to the diffusion-controlled limit, extension L d = 30 000 1 , 400 2 , 50 3 , 25 4 , tia ( 0 ) = tib(0)) 


d 

1 

2 

3 

4 

0S 

0.516011 



0.98217 


0.4966i2 

0.9686a 

0.9390a 

1.0028 9 

poo 


0.9504a 

0.9736a 

0.9575a 

1.0025 9 


Table 2: Numerical results of A + A —> 0, see figure (|7). As in the previous table the subscript 
s denotes the order of magnitude of the reaction constant. Note, that in d, = d c = 2 logarithmic 
corrections have been performed. 


causes an additional overhead by a factor of 5-10. 

The typical speed of a common Ultrix or Alpha 80 Me RISC workstations leads to some 
10 4 — 10 5 Markoff-operations per second. For some 10 7 annihilation processes we need some 
10 minutes to a few hours. 

We claim our algorithm to be superior to any lattice gas simulation, because we are not 
restricted to one molecule per cell. For k —> 00 the lattice gas is included as limit, for k —> 0 
the simulation degenerates to the homogeneous case. An upper limit for the size ui of the 
reaction cells is given by the comparison of reaction and diffusion-time-scales. 

The computing time for more complex RDS will increase linearly with the number of 
reactions. In our simple cases the memory overhead also is about a factor of 2, and can be 
compared to two additional components of the RDS. However, after we have implemented 
an event handler which is completely independent of physical problem and the underlying 
topology, further simulations of RDS are realized with little effort. The event handler can be 
easily integrated into any program, because the only necessary adaptions are the selection of 
the reactions in a single cell. 

6.1 Future extensions 

Due to the generality of the method, we see a lot of possible applications. 

6.1.1 Obvious generalisations 

Further extensions are possible. If its neighbourhood influences the reactions of a cell, i.e. if 
we want to describe surface dynamics or the fluctuations of a field, the general method keeps 
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unchanged. Only a slight modification of the update of the neighbour cells must be taken into 
account. 

The simplest assumption, the chemical reactions being independent of ft, may be com¬ 
pleted by reactions depending on the single cell or a cell and its topological neighbourhood. 
With regard to these conditions we only have to supplement further indices if necessary. The 
transition probabilities, which do not only depend on the single cell n, but also on its neigh¬ 
bourhood may be denoted by the symbol This way we can indicate the possibility 

that a step depends on the contents of more than one single cell. For those types of events the 
symbol M is in contrast to A, which shall be left to local rates. Therefore depends 

on the cells n,m,m'.. whereas Ah is restricted to n. If we provide Ah^^, with additional 
indices, this symbol serves as a denotation for a transition changing the contents of several 
cells n, m, m' with a rate depending only on the single cell n. This is the manner we have 
described diffusion in Gillespie’s algorithm by a reaction changing the contents of two cells, 
whose rate A a fim only depends on n alone. 

To handle a reaction as a Markoff-event by the event handler, it is always assigned to its 
first cell specified. Thus AT and both are assigned to cell n. 

Q ft = A ^ + Mfi 

= Ah A +• • •++• • •. (58) 

i i i i i 

The difference of the dependencies on one or more cells is unimportant for the principal 
selection of a reaction step. Again an arbitrary event is unambiguously denoted by i and n 
and its rate Q\ ; . More precisely, this symbol now unifies the description of Ah, Ah^, . 

For a transition touching further cells m, rh! the bookkeeping had to be done for these cells, 
too. 

For RDS the case of a reaction depending on several cells does not attend, we have presented 
this generalised formalism to show, that even more complex dynamics could be handled. This 
way it is possible to introduce events depending on gradients into the framework of chemical 
reactions. This opens an outlook not only to the simulation of surface dynamics, but points 
out to new approach to field theoretical models on computers. 

We want to proceed to simulate more complex behaviour, with a variety of possible chemical 
reactions. 

6.1.2 Automatically generated code 

Therefore, we are writing a code generator, generating C source code by a script to automatise 
the process from the problem definition to the program running. As we have seen, more 
complex reactions are difficult to implement because for any practical purpose, any efficient 
code expands from one line definition to one page C source code. This will be the subject of 
a future paper. 

We expect even much more interesting stochastic effects for more complex systems or in 
more complicated topologies. Simulations become essential, because by the lack of a simple 
scaling ansatz we may lose control over the fluctuations in a complex system, especially if the 
numbers of molecules in an individual cell are small. 
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: Scaling approach for A + A —> 0 In the diffusion-controlled limit, the mean field behaviour 
cannot correctly describe the role of fluctuations. The approach of the scaling theory is to 
derive arguments connecting fluctuations to a length scale £. The length £ is related to a 
typical time scale t £ by standard diffusion-scaling = £ 2 /D. 

Kang and Redner have argued [!(J, that the initial mean distance £ = a aw of a A molecule 
to its next neighbour is proportional to n(0) _1 ^ d . Therefore, the time-scale tnn to pass this 
distance is scaled by the standard diffusion-scaling of space and time according to 


with the asymptotic properties 


tnn — 

a NN 

D 

= n{0)~i. 

(59) 

n as g(x) 



n(t) = 

n(0) • g 

(—)• 

(60) 



\tnn / 


x < 1 

9(x) 

= 1, 

(61) 

x » 1 

9(x) 

oc x^ 13 . 

(62) 
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To compute the scaling exponent (3, we need the additional assumption, that the long time 
behaviour is independent of the initial concentration, thus we obtain 

(Dt)^n(t) = constant. (63) 

In d c the scaling and the mean-field exponents meet, thus 

d c = 2. (64) 

: Scaling approach for A + B —> 0 The dynamics of this system are governed by domains, 
which are defined as local areas where either the A or the B species dominates. A typical 
picture obtained by simulations in d = 2 is shown in Fig. 

Because A and B annihilate as pairs, the difference 


An = riA(t) -ns(f) = ua( 0) — ub( 0) = constant (65) 

is conserved. In a volume of linear extension £ the initial number of molecules of each 
species may vary due to fluctuations at time t = 0 

N a (0) = {N a ) ± y/{Nk) = n Q (0)C d ± y/n a { 0)( d . (66) 


Some species may locally be in majority even it is global in minority i.e. Na > Nb even if 
riA < riB- Thus by fluctuations local inversions are possible. The largest possible volume £ d , 
in which an inversion may occur, is estimated by equating the numbers of A and B molecules 


Inserting 


0 = Na (0) — Nb (0) = | VNa + V Nb | • 

]) we solve © obtaining the length scale 


e = 


\J ? ia ( 0 ) - \Jn B ( 0 ) 


(67) 


( 68 ) 


By diffusion-scaling £ is related to the time required to cross a domain 

k = y r (69) 

assuming both diffusion-constants being equal, i.e. D = Da = Db■ The scaling ansatz 

n a (t) = Cat - ' 1 f a ( 70 ) 

requires the existence of scaling functions f a {x), whose exponential decay dominates the alge¬ 
braic decay t 1 for t > t%. For short times the scaling functions are assumed to be constant, 
thus 


Inserting f a into (^) we get 


x < 1 
x » 1 


U(x) 

fix) 


= 1 , 

oc exp(— cx). 


Ca ' fA 



-°sh {£) 


i \J tia(0) - yj ns(0) • i t 1 . 


(71) 

(72) 


(73) 


The left side of the equation only depends on the ratio t/t^, therefore this statement yields 
for the right side of this equation, too, thus leading to the algebraic exponent 



(74) 


Although in the limit ua( 0) —> ua{ 0) the factor 


\/™ a ( 0 ) - yJn B ( 0 ) 


0 in equation (n3) 


disappears, the argumentation is not affected and the scaling keeps valid 10- 
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